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We consider a bilayer of dipolar bosons in which the polarization of dipoles are perpendicular to the 
planes. Using accurate static structure factor S(q) data from hypernetted-chain calculation for single layer 
dipolar bosons we construct effective screened interactions for intralayer particles. We adopt the random- 
phase approximation for interlayer interactions. We study the instability of the homogeneous bilayer system 
against the formation of density waves by investigating the poles of the density-density response function. 
The dispersion of collective modes of this system also signals the density-wave instability. We also inves¬ 
tigate the effect of counterflow on the collective mode dispersion and on the density-wave instability and 
discuss the dissipationless superfluid drag effect in the presence of a background velocity. 
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I. INTRODUCTION 

Experimental progresses in trapping and cooling atoms 
with large magnetic moments and polar molecules, opened 
up a new and interesting area of exploring quantum many- 
body systems with large and anisotropic dipole-dipole in¬ 
teractions [UHl- Bose-Einstein condensation (BEG) of 
polar molecules [RHtH. and atoms with large permanent 
magnetic moments QlHTl has been observed. Dipolar 
fermionic gases have been cooled down near to their 
ground-state as well 11111 1 1211 . 

In bulk geometries, the attractive part of the dipole- 
dipole interaction could in principle lead to instabilities, 
like BEG collapse II13II or chemical reactions between parti¬ 
cles []3l]. Therefore, it is usually favorable to confine the 
dipolar gases into quasi-two or one-dimensional geome¬ 
tries, and use an external electric or magnetic field (de¬ 
pending on the nature of dipoles) to polarize all dipoles in 
the same direction. Layered structures are another config¬ 
uration of great interest which one can tune the attractive 
interactions and pairing between different layers without 
the fear of having chemical reactions. 

While the stripe or density-wave phase is naturally ex¬ 
pected in an isolated 2D system of tilted dipolar bosons 11141 
EH and fermions GMH due to the anisotropy of 
the dipole-dipole interaction, this instability has been 
the subject of much dispute in the limit of perpen¬ 
dicular dipoles, where the inter-particle interaction is 
isotropic. While mean-field approximation II16I1 as well 
as density-functional theory (DET) II19II and Singwi-Tosi- 
Land-Sjdlander (STLS) lIlSll calculations all predict stripe 
phase formation at relatively low interaction strength 
for two-dimensional dipolar fermions, quantum Monte 
Garlo (QMG) simulations find that the stripe phase never 
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becomes energetically favorable, up to the liquid-to- 
solid phase transition for perpendicular bosons lllSlI and 
fermions ||20|| . 

In double-layer structures, the fermionic systems at¬ 
tracted enormous attention. Eor example the ground state 
properties and instabilities have been studied within the 
Hartree-Eock II21H23II as well as STLS methods II24II . The 
QMG simulation has been employed to study the crossover 
from BEG to Bardeen-Gooper-Schrieffer (BGS) state II25II 
too. Bosonic bilayers, on the other hand, received slightly 
less attention. Hufnagl and Zillich [l26ll used h5q3ernetted- 
chain approximation (HNG) to calculate the ground-state 
quantities of a bilayer system of tilted dipolar bosons. Then 
using the correlated-basis function method they obtained 
the dynamical properties. More recently, the competition 
between single-dipole and dimer condensation in a bilayer 
of perpendicular dipolar bosons have been investigated us¬ 
ing QMG method by Macia et al. II27II They have observed 
that at strong interlayer coupling, the pair superfluidity 
dominates over the single-particle superfluidity. 

In this work we study two parallel layers of dipolar 
bosons, whose moments are aligned perpendicular to the 
2D-plane. Therefore, both the intralayer and interlayer in¬ 
teractions are isotropic. The bare intralayer interaction is 
purely repulsive while the bare interlayer interaction could 
be either repulsive or attractive, depending on the in-plane 
separation of two dipoles [see, Eqs.Q and below]. We 
investigate the possibility of instability of a homogeneous 
liquid towards the formation of inhomogeneous densities, 
i.e. density waves. For this purpose we look at the poles 
of the static density-density response function. The effec¬ 
tive intralayer interactions are obtained from an accurate 
HNG results for the static structure factor of an isolated 
2D layer II28II . combined with the fluctuation-dissipation 
theorem II29II . and the interlayer interactions are treated 
within the random-phase approximation (RPA) Il30ll . We 
also find the full dispersion of in-phase and out-of-phase 
collective modes (i.e. zero-sound modes) from the poles of 
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the dynamical density-density response function. We have 
also included a finite counterflow between two layers and 
studied its effects on the density-wave instability and col¬ 
lective modes. We finally address the dissipationless super¬ 
fluid drag in the presence of a finite background flow. 

A similar study of the instability of a homogeneous liquid 
with respect to the inhomogeneous phase of charge density 
wave has been studied in a variety of quantum charged 
systems ranging from single-layer electron gas II31II to 
electron-electron and electron-hole double-layers II321135II 
to charged Bose systems II36II and superlattices II37II . 

The rest of this paper is organized as follows. In sec¬ 
tion]^ we introduce the density-density response function 
of our system and explain how effective intralayer interac¬ 
tion could be obtained from the static structure factor. In 
section]^ we calculate the density-wave instability and the 
dispersion of the collective modes. Finally in section 
summarize and conclude our main findings. 


IV 


we 


II. DENSITY-DENSITY RESPONSE FUNCTION AND THE 
EFFECTIVE INTERACTIONS 


We consider two identical two-dimensional planes of 
bosonic dipoles, separated by distance d. All dipoles are 
assumed to be polarized perpendicular to the planes. The 
bare intralayer and interlayer interactions respectively read 


and 


f^dd f 
V(r)=— — 

s'- A S’ 

471 


Cdd r^-2d^ 

471 (r^ + ’ 


(1) 


( 2 ) 


where Qd is the dipole-dipole coupling constant and r is 
the in-plane distance between two dipoles. After Fourier 
transformation one finds ||22|| 


V.(q)= — 


and 


3V2^w 


, (3) 




(4) 


Here erfc(x) is the complementary error function and w is 
the short distance cut off introduced to heal the divergence 
of Fourier transform of the intralayer interaction. 

In this work we are interested in the density-wave in¬ 
stabilities and collective density modes of this double layer 
structure. For this we begin with the following matrix equa¬ 
tion for the density fluctuations II30II 


j 


where 5n^[q,co) is the density fluctuation in layer i (i = 
1,2). Vl“‘(q,£u) is the external potential applied to layer 


j and Xij(q,co') is the density-density response function, 
whose matrix form reads 

X~^[q,co) = n~\q,co)-W^^\q,co) . ( 6 ) 

Here Ylij[q,co) = 5yn;(q,m) is the non-interacting 
density-density response function of layer i, and W^®(q, co) 
is the dynamical effective potential. For symmetric bilay¬ 
ers we have Yl^[q,co) = n(q,£u) (same for both layers), 
and W^®(q, m) = 5;^-W,(q, m) -F (1 - S^OW^Cq, co), where 
Ws(q,a)) [Wjj(q,(u)] is the effective interaction between 
dipoles in the same [different] layers. 

Eigenvalues of the density-density response matrix 
2[q,co) are 




n(q,m) 

l-n(q,m)W±(q,m) ’ 


(7) 


where W±(q, co) = WjCq, co) ± W^[q,co) are the symmetric 
and antisymmetric components of the effective potentials. 

The non-interacting density-density response function 
n(q,£u) of a two dimensional system of bosons is analyt¬ 
ically known at zero temperature 


n(q,m) 


2nEq 

{hco + i0+)^ - 


( 8 ) 


where n is the particle density in one layer and = 
h2q^/(2m) is the single-particle energy. The exact form 
of the effective potentials are not known, and one has to 
resort to some approximations. In the celebrated random 
phase approximation II30II . the effective potentials are re¬ 
placed with their bare values. In this work we will use the 
bare interlayer potential V[j(q) from Eq. On the other 
hand, the bare intralayer potential, in q-space, ([^ has an 
artificial cut-off dependence. In order to overcome this 
problem, we will use the fluctuation-dissipation theorem 
to find an approximate expression for the interlayer poten¬ 
tial II29II . At zero temperature the fluctuation-dissipation 
theorem reads II30II 


S(q) 


Htc f°° f n[q, co) 

— dm 3m -T;-r 

^ Jo l'i--W,[q,co)niq,co) 


(9) 


Now, neglecting the frequency dependence of the effective 
potential i.e. replacing W^Cq, co) with a static and real func¬ 
tion Ws(q), we can perform the frequency integral in Eq. ^ 
analytically and obtain 




2n 



( 10 ) 


Here, S(q) is the static structure factor of an isolated 2D 
dipolar Bose liquid, which can be obtained very accu¬ 
rately e.g., from QMC simulations 11381139II or HNC cal¬ 
culations Il28ll . The effects of exchange and correlation, 
which are entirely ignored in the RPA, are partly included 
in Eq. (101 through the static structure factor. We set the 
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interlayer part of the effective interaction to the bare inter¬ 
layer interaction i.e., W^(q) = V^{q) since a QMC or HNC 
calculation for interlayer static structure factor is not avail¬ 
able. Such an approximation is equivalent to RPA and we 
surmise it will be adequate for large enough layer separa¬ 
tions. 


III. DENSITY-WAVE INSTABILITIES AND COLLECTIVE 
MODES 


Density-wave instabilities and collective modes could be 
obtained from the poles of the density-density response 
function given in Eq. Q, respectively in the static and dy¬ 
namical limits, or equivalently from the solution of 

l-n(q,m)W±(q,m) = 0, (11) 


Density-wave instability: In the static limit, the non¬ 
interacting density-density response of Eq. reduces to 

2n 

n(q) =-, (12) 

which together with Eq. ( [Tl] ), gives 
2n 

1+—W±(q) = 0. (13) 

Now, using the bare interlayer potential Q and the effec¬ 
tive intralayer potential of ( [Tol ) in Eq. ( [T^ we find 

q = ±8nnrQS^[q)e~^‘^ . 


Here Vq = mQd/(47rfi^) is a characteristic length scale. 
As the static structure factor is positive by definition, the 
above expression with minus sign will not have any solu¬ 
tion which means that no density-wave singularity in the 
out-of-phase channel (i.e., X- = 0) appears. On the other 
hand, in the in-phase channel (i.e., X-^- = 0) one can find 
instabilities for suitable values of the interaction strength 
and layer spacing from the solutions of Eq. ( [T^ w ith the 
positive sign. Numerical investigation of Eq. (141 reveals 
that (see, Fig.[^ at every value of the coupling constant 
Y = nr^, there is a critical layer spacing d^, below which 
one can find density-wave instability at a finite value of the 
wave vector q^. 

We note that for an isolated single layer, one has W^[q) = 
0, and the criteria for the density-wave instability becomes 


0 = 1- W,(q)n(q) = , (14) 

S^(q) 

which has no solution at any finite q. Therefore, within 
the approximations we use here, no density-wave instabil¬ 
ity is expected to happen in an isolated two-dimensional 
system with purely repulsive dipolar interaction. In agree¬ 
ment with the QMC findings lIlSll . 

The behavior of static density-density response functions 


z±(q) 


4mn 1 


(15) 



FIG. 1. The critical layer separation versus coupling constant y 
for bilayer dipolar bosons. The purple region shows the homoge¬ 
neous superfluid (SF) phase and the green one is the region with 
density-wave (DW) instability. 


are also illustrated in Figs.[^and[^ As the layer separation 
is lowered to the critical distance, a singularity in 
emerges (see, Fig.|^, but the antisymmetric component of 
the density-density response function remains finite. 

Collective modes: In order to find the collective modes, 
we study the singularities of the dynamical response func¬ 
tions ;y±(q,a)). If we approximate the dynamical effective 
interactions with static and real functions, we find from 
Eq. 

hfcoliq) = £q+ 2neqW±(q) . (16) 

Again, replacing the effective interlayer potential W^[q) 
with the bare interaction V'd(q), and the effective intralayer 
potential W^q) from Eq. ( [T^ , the dispersion of collective 
modes read 








TnC^^qe 


(17) 


Note that the first term on the right-hand-side of this equa¬ 
tion is the Bijl-Feynman excitation spectrum of a single 
layer dipolar Bose liquid II28II . In the long wavelength limit, 
as the static structure factor vanishes linearly [S(q ^ 0) oc 
q], we find 


<w±(q) v,q -h 0[q^) , (18) 

where is the sound velocity. Note that unlike the bilayer 
charged boson system II40II , both collective modes of a bi¬ 
layer system of dipolar bosons have acoustic nature and 
are degenerate at long wavelengths. Using the numerical 
results for the static structure factor from Il28ll in Eq. ([T3i, 
the full dispersion of the collective modes could be readily 
obtained. The results for a)±(q) and single-layer collective 
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FIG. 3. The antisymmetric component of the static density-density 
response function ;f_(q) as a function of the dimensionless wave 
vector q/v^ for several values of the layer spacing d at 7 = 16. 


Counterflow: Now we consider background velocities of 
Vj and V 2 in the first and second layers respectively. The 
collective modes and density-wave instabilities in the pres¬ 
ence of these background velocities could be easily ob¬ 
tained from the poles of the total density-density response, 
after replacing co with m — V; • q in the noninteracting den¬ 
sity response of layer i, n;(q,m). The background flows 
could be decoupled into the center-of-mass V = (v^ -l-V2)/2 
and counterflow v = (v^ —wf)!! components. As the ef¬ 
fect of the center-of-mass flow could be simply understood 
in terms of a Galilean boost, we focus on the counterflow 
part, and look for the solutions of the following equation 


FIG. 2. The symmetric component of the static density-density 
response function /+(q) as a function of the dimensionless wave 
vector ql^/n for several values of the layer spacing d. As d ap¬ 
proaches the critical spacing d„ a singularity at finite q emerges 
in the density-density response function. Top and bottom panels 
are for two different values of the coupling constant, 7 = 16 and 
7 = 64, respectively. 


mode cUjjCq) = Eg/[hS[q)] [in units of Eq = hfn/{2m)] are 
presented in Fig.j^for a fixed value of the coupling constant 
7 and for different values of the layer separation d. We find 
that the critical layer spacing is = 0.484/-/n for 7 = 16, 
at which point the symmetric mode £u+ touches zero and 
becomes soft. This occurs at the same q-value that the sym¬ 
metric density-density response function x+ diverges (c.f. 
Fig.0. 

In Fig.j^the full dispersions of collective modes are pre¬ 
sented at different values of the coupling constant 7 . We 
kept the layer spacing at the critical value of the layer 
dj. = 0.5507/^/n for 7 = 32 and observed a similar mode 
softening with increasing 7 . 


n(q,m-vq) 

-Wd(q,m) 


Ws(q,m) 


-Wd(q,m) 


n(q,m + v-q) 


- W.(q,m) 


0 . 


(19) 


In the static limit (m = 0), we find 

2 n f hv-q\^ 

1 + — W±(q) - -^ = 0 . ( 20 ) 

V ; 


Now, using the same approximation for the intra- and inter¬ 
layer potentials as before, we arrive at 


q = ±8nnrQS^[q)e + 


4m^v^ cos^ 6 
h^q 


S^(q) , 


( 21 ) 


where 9 is the angle between counterflow direction and 
the q vector. Eq. (211 clearly shows that at a finite counter- 
flow, the density-wave instability in the symmetric channel 
is facilitated in the directions parallel to the flow. 

Similarly, for the collective modes in the presence of a 
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q/Vn 

FIG. 4. Dispersion of symmetric and antisymmetric co_ modes 
at a fixed value of the coupling constant y = 16, and for differ¬ 
ent values of the layer separation: d = 0.49 (top), d = 0.484 
(middle), and d = 0.47 (bottom). The dashed line represents 
single-layer collective mode cUj/Cq). Note that d = 0.484 is the 
critical value of the layer separation of the formation of density- 
wave instability for 7 = 16 . 



q/v4i 


FIG. 5. Dispersion of collective modes of a bilayer system of dipo¬ 
lar bosons for different values of the coupling constant 7 = 16 
(top), 7 = 32 (middle), and 7 = 64 (bottom) and for a fixed 
value of the layer spacing d = 0.5507. The dashed line represents 
single-layer collective mode cUjjCq). 
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main effect of lowering the symmetric mode and 

enhancing the antisymmetric mode co_(q). As v increases 
symmetric mode becomes zero at two values of q indicating 
instabilities according to Eq. ( [2T] ). 

Zero-point energy and drag effect: With the full dispersion 
relations, we are able to find the change in the zero-point 
energy, due to the finite counterflow ||41ll 

Xi ['^aCq.v)-m„(q,0)] , (26) 

q,a=± 


in which the difference between collective modes with and 
without counterflow are summed. Now, using the leading 


order results in v, from (241 and assuming that the static 
structure factor is power expandable at small q, e.g.. 


S(q - 0) = 


hq 

2 mv, 


+ PiQ. + ftq + ••• > 


(27) 


FIG. 6. Dispersion of collective modes of a bilayer system of dipo¬ 
lar bosons at finite counterflow velocity v and several values of 
the layer separation d. The interaction strength is y = 16. Note 
that we specialize to the case 6 = 0 in Eq. ( [23[ l, i.e., when the 
wave vector q is parallel to v. 


counterflow, we obtain 


h^m^(q,v) = e; 


2n 

1 + — W/q) 

£n 


-h [hv-qY 


/ 

2n 

1 n^W^Cq]-h (Hv qY 

1 + —W,(q) 


( 22 ) 


We note that the dispersions of collective modes become 
anisotropic in the presence of a finite counterflow. Finally, 
replacing the effective interaction WjCq) from Eq. ( [T^ , and 
the bare interlayer interaction V^Cq) from Eq. (Q, we find 


h^m^(q,v) = 


_s(q)_ 


-b [hv-qY 


=Fe„Wn2C?,q2e -p 


2hv •q 1 ■ 
_ S(q) 


(23) 


To leading order in the counterflow velocity v, we find for 
the dispersions 


we find 





hC^^n^v^ cos^[6) 
32m^v^ 

5 


q3 + 0[qY 


(28) 


Expansion of each dispersion alone depends on the higher 
order coefficients of S(q), but interestingly they cancel out 
in the summation of two modes (up to q^ terms). Integrat¬ 
ing this expression up to q ~ 1/d, we obtain 


AEzp 


6407rm^Vj^d^ 


(29) 


We note that the zero-point energy is negative, meaning 
that a finite counterflow lowers the energy of the system, 
depends on the square of the velocity and scales as d~^ 
with the separation of two layers. It is interesting to note 
that zero-point energy for dipolar fermions yields d~^ de¬ 
pendence on the separation distance. II42II Previous cal¬ 
culations II40II on charged systems interacting via the 1/r 
Coulomb potential have shown that the bilayer separation 
dependence of the zero-point energy is different for elec¬ 
trons and charged bosons too. In Fig.j^we show the nu¬ 
merical results of AE^p as a function of the counterflow 
velocity v at d = 1 / -/n for various values of the coupling 
strength j. 

We now construct the free energy E, by adding the ki¬ 
netic energies of the fermions in each layer 


m±(q,v)f%J m±(q,0)-|- 


mS(q)v^cos^(0) 

h^lTa(q) 


where 



2 ■ 
a(q) 




(24) 


a(q) 


SnnrQ 

q 


S\q)e-’i‘^ . 


(25) 


In Fig.j^we illustrate the effect of counterflow velocity 
V on the collective modes for a bilayer dipolar fluid at y = 
16 and d = l/\/n. We note that the counterflow has the 


E = -nm(Vj2-|-v|)- 


hC^^n 


2560nm^\ 




(30) 


where we have reverted to use the individual velocities in 
each layer. We find the current densities in layer 1 and layer 
2 , calculated from = dF /d V; (i = 1,2), to be 


Ji = 


12807Tm^vf 


Vi + 


12807rm^v,^d^ ^ ’ 


J2 — 


hCdd^ri^ 


1280nm^vfd^ 


Vi + nm 


hc^dV A 

12807rm^v^^d^ J ^ * 


( 31 ) 
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FIG. 7. Change in the zero point energy AE^p as a function of 
counterflow velocity v and for several values of the interaction 
strength y. The layer separation is d = \l sfn. 

The expressions given in Eq. ( [31] ) demonstrate that the su¬ 
perflow in the first (second) layer depends on the super¬ 
fluid velocity on the same layer as well as that of the sec¬ 
ond (first) layer. This is the dissipationless superfluid drag 
effect well known in two-component superfluids ||40tf45|| 
which has been discussed for a variety of related systems. 

IV SUMMARY AND CONCLUSIONS 

In summary, we have investigated the instability of a ho¬ 
mogenous bilayer system of perpendicular dipolar bosons 


towards density waves. An accurate HNC results for the 
intralayer static structure factor is used together with the 
fluctuation-dissipation theorem to extract a static intralayer 
effective potential and the random phase approximation is 
employed for the interlayer interaction. We have observed 
that at any intralayer coupling strength, there is a critical 
layer spacing below which the homogenous in-plane den¬ 
sity becomes unstable. The full dispersion of the in-phase 
and out-of-phase zero-sound modes of the bilayer system 
have been calculated too. We observed that two modes be¬ 
come degenerate in the long wave-length limit. 

Furthermore, we have studied the effects of a finite 
counterflow on the density-wave instability and collective 
modes. Counterflow lowers the zero-point energy and 
drive the homogenous system towards the formation of 
density waves. We have also studied the dissipationless 
drag between the two layers of dipoles similar to other bi¬ 
layer systems. 

Finally, we should note that in the limit of closely sepa¬ 
rated layers, improvements beyond the RPA in the effec¬ 
tive interlayer potential, like the inclusion of exchange- 
correlation effects or the binding of two dipoles from dif¬ 
ferent layers might be necessary. Dynamical effects and 
frequency dependance of the effective potentials would be¬ 
come important in the strongly-correlated regime too. 
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